En este trabajo realizaremos un intento de predecir las hipotecas de la provincia de Cáceres durante en el año 2019 empleando un modelo ARIMA en el contexto de las series temporales.

Carga de librerías y datos

Importamos todas las librerías y asignamos el directorio de trabajo, además de cargar el fichero csv que contiene nuestros datos.

library(ggplot2)
library(plotly)
library(tseries)
library(MASS)
library(forecast)
library(lmtest)
library(caschrono)
library(lubridate)
library(timeDate)
library(tsoutliers)
library(dplyr)

# Se limpia el espacio de trabajo:

rm(list = ls())

# Se asigna el directorio donde están los datos

setwd("C:/Users/Javier/Documents/masterAFI/16. Analisis de series temporales/01. Series temporales/javier_herraez_albarran_series_temporales/")
# setwd("C:/Users/jherraez/Documents/masterAFI/16. Analisis de series temporales/01. Series temporales/practica/")

datos <- read.csv("_11_CACERES.csv")

Dividimos en train y en test, y convertimos los datos en un objeto de tipo time series para poder aplicar sobre ellos las funciones de R.

datos$FECHA = as.Date(datos$FECHA,format='%d/%m/%Y')
datos.train <- subset(datos, FECHA<=as.Date('01/12/2018',format='%d/%m/%Y'))
datos.test <- subset(datos, FECHA>as.Date('01/12/2018',format='%d/%m/%Y'))

datos.train.ts <- as.ts(datos.train$TARGET, frequency = 12)
datos.test.ts <- as.ts(datos.test$TARGET, frequency=12)

Representación gráfica de los datos

Elaboramos un gráfico que nos permita ver la evolución de las hipotecas en los últimos años.

names(datos.train)
## [1] "FECHA"  "TARGET"
graficoInicial <- ggplot(aes(x= FECHA, y = TARGET), data = datos.train) +
  geom_line(color = '#d84519', size = 1) + 
  xlab('FECHA') + ylab('Matriculaciones')

ggplotly(graficoInicial)

Vemos como se han reducido considerablemente las hipotecas en estos años más próximos, aunque se puede seguir intuyendo algo de periodicidad.

Estacionalidad

Estacionalidad en varianza

Evaluamos si es necesario transformar la serie para hacerla estacionaria en varianza. Para ello realizamos un boxcox.

box_cox <- boxcox(TARGET ~ FECHA,
                  data = datos.train,
                  lambda = c(0, 0.5, 1))

lambda <- box_cox$x[which.max(box_cox$y)]
lambda
## [1] 0.1616162

Como lambda toma un valor cercano a cero, tomamos logaritmos.

datos.train$log_target=log(datos.train$TARGET)
datos.test$log_target=log(datos.test$TARGET)
datos.train.ts <- as.ts(datos.train$log_target, frequency = 12)
datos.test.ts <- as.ts(datos.test$log_target, frequency=12)

Estacionalidad en media

Probamos a hacer el test de Dickey-Fuller aunque no confiamos demasiado en su resultado, ya que es poco exigente.

adf.test(datos.train.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  datos.train.ts
## Dickey-Fuller = -1.6659, Lag order = 5, p-value = 0.7159
## alternative hypothesis: stationary

Todo parece indicar que va a ser necesaria una diferencia regular, pero nosotros ajustamos directamente sin nada.

Ajuste de un modelo SARIMA a la serie

Observamos los gráficos f.a.s. y la f.a.p. para realizar el primer ajuste.

acf(datos.train.ts, lag.max = 48, xlab = "Retardo",
    main= "Función de autocorrelación simple")

pacf(datos.train.ts, lag.max = 48, xlab = "Retardo",
     main = "Función de autocorrelación parcial")

Lo habitual es empezar con una estructura AR(1). La fuerte correlación de orden 1 en la f.a.p parece reforzar esta decisión, por lo tanto, se propone un ARIMA(1,0,0). Vemos los coeficientes de las variables que obtenemos al realizar este primer modelo.

ajuste1 <- Arima(datos.train.ts,
                 order = c(1,0,0),
                 seasonal = list(order = c(0,0,0), period = 12),
                 method = "ML")
ajuste1
## Series: datos.train.ts 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.9059  5.5001
## s.e.  0.0303  0.2199
## 
## sigma^2 = 0.09146:  log likelihood = -42.67
## AIC=91.34   AICc=91.47   BIC=101.11
coeftest(ajuste1)
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ar1       0.905855   0.030348  29.849 < 2.2e-16 ***
## intercept 5.500117   0.219941  25.007 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Probamos también a ver las correlaciones de estas.

cor.arma(ajuste1)
##                   ar1   intercept
## ar1        1.00000000 -0.02691084
## intercept -0.02691084  1.00000000

Vemos que en valor absoluto estarían por debajo de 0.8 por lo que no las consideremos correladas.

Mediante el test de LJung-Box estudiamos si hay ruido blanco.

Box.test.2(residuals(ajuste1),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard  p-value
## [1,]      6 6.49e-06
## [2,]     12 5.80e-07
## [3,]     18 2.54e-06
## [4,]     24 1.10e-07
## [5,]     30 6.80e-07
## [6,]     36 8.00e-08
## [7,]     42 4.30e-07
## [8,]     48 3.00e-07

Buscamos obtener unos p-valores mayores a 0.05, por lo que como vemos no tenemos ruido blanco.

Seguimos ajustando el modelo, así que observamos los gráficos f.a.s. y la f.a.p. para realizar un segundo ajuste.

acf(ajuste1$residuals, lag.max = 48, xlab = "Retardo",
    main= "Función de autocorrelación simple")

pacf(ajuste1$residuals, lag.max = 48, xlab = "Retardo",
     main = "Función de autocorrelación parcial")

Según lo que observamos en las gráficas, como destaca la línea del 1 en la gráfica de autorrelación simple, ajustamos un modelo ARIMA(1,0,1).

ajuste2 <- Arima(datos.train.ts,
                 order = c(1,0,1),
                 seasonal = list(order = c(0,0,0), period = 12),
                 method = "ML")
ajuste2
## Series: datos.train.ts 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    mean
##       0.9947  -0.7254  5.5238
## s.e.  0.0054   0.0434  0.5351
## 
## sigma^2 = 0.06019:  log likelihood = -2.52
## AIC=13.04   AICc=13.26   BIC=26.07
coeftest(ajuste2)
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.9947456  0.0054474 182.610 < 2.2e-16 ***
## ma1       -0.7254036  0.0434252 -16.705 < 2.2e-16 ***
## intercept  5.5238319  0.5351063  10.323 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vemos que el coeficiente del AR1 + 2 * Std. Error pasaría del valor 1, así que es claro que necesitamos diferenciarlo. Pasamos entonces a un modelo ARIMA(0,1,1).

ajuste3 <- Arima(datos.train.ts,
                 order = c(0,1,1),
                 seasonal = list(order = c(0,0,0), period = 12),
                 method = "ML")
ajuste3
## Series: datos.train.ts 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.7276
## s.e.   0.0426
## 
## sigma^2 = 0.05974:  log likelihood = -1.79
## AIC=7.59   AICc=7.65   BIC=14.09
coeftest(ajuste3)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.727557   0.042638 -17.064 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hacemos el test de Ljung-Box

Box.test.2(residuals(ajuste3),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard    p-value
## [1,]      6 0.67336896
## [2,]     12 0.23278469
## [3,]     18 0.25043194
## [4,]     24 0.04279890
## [5,]     30 0.09096946
## [6,]     36 0.03682012
## [7,]     42 0.04956558
## [8,]     48 0.07392056

Vemos que ya casi todos los valores están por encima del 0.05, aún así seguimos ajustando.

acf(ajuste3$residuals, lag.max = 48, xlab = "Retardo",
    main= "Función de autocorrelación simple")

pacf(ajuste3$residuals, lag.max = 48, xlab = "Retardo",
     main = "Función de autocorrelación parcial")

Según las gráficas, vemos que podemos proponer un SARIMA(0,1,1)x(1,0,0)[12]. Construyamos este modelo..

ajuste4 <- Arima(datos.train.ts,
                 order = c(0,1,1),
                 seasonal = list(order = c(1,0,0), period = 12),
                 method = "ML")
ajuste4
## Series: datos.train.ts 
## ARIMA(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.7430  0.2382
## s.e.   0.0445  0.0729
## 
## sigma^2 = 0.05668:  log likelihood = 3.35
## AIC=-0.71   AICc=-0.58   BIC=9.05
coeftest(ajuste4)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ma1  -0.743005   0.044523 -16.6881 < 2.2e-16 ***
## sar1  0.238164   0.072912   3.2664  0.001089 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Medimos también correlaciones.

cor.arma(ajuste4)
##             ma1       sar1
## ma1   1.0000000 -0.1216682
## sar1 -0.1216682  1.0000000

Y volvemos a hacer el test de Ljung-Box para ver si ya obtenemos ruido blanco.

Box.test.2(residuals(ajuste4),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard   p-value
## [1,]      6 0.7123801
## [2,]     12 0.9055981
## [3,]     18 0.8138609
## [4,]     24 0.5663195
## [5,]     30 0.7499105
## [6,]     36 0.6405039
## [7,]     42 0.6366659
## [8,]     48 0.7534805

Vemos que tenemos valores por encima de 0.05 para todos, así que parece que podemos decir que tenemos ruido blanco.

Análisis de intervenciones

Crearemos una función para crear variables explicativas en series mensuales. Esta función vista en clase añade alguna variante para añadir festivos propios de Cáceres.

calculoExplicativasCalendarioCaceres <- function(variableFecha, domingoYFestivosJuntos){
  
  #######################################
  #     Creación de todas las fechas    #
  #######################################
  
  if (month(max(variableFecha)) %in% c(1,3,5,7,8,10,12)) {
    diasHastaFinMes <- 30
  } else if (month(max(variableFecha)) %in% c(4,6,9,11)) {
    diasHastaFinMes <- 29
  } else if (year(max(variableFecha))%%4==0) {
    diasHastaFinMes <- 28
  } else {diasHastaFinMes <- 27}
  
  todasLasFechas <- data.frame(fechas=seq(min(variableFecha),
                                          max(variableFecha)+diasHastaFinMes,
                                          by="days"))
  
  #######################################
  #     Cálculo de la Semana Santa      #
  #######################################
  
  domingoResurrecion <- as.Date(Easter(year(min(variableFecha)):year(max(variableFecha))))
  lunesPascua <- domingoResurrecion + 1
  sabadoSanto <- domingoResurrecion - 1
  viernesSanto <- domingoResurrecion - 2
  juevesSanto <- domingoResurrecion - 3
  
  # Se unen y ordenan todos los días que forman la Semana Santa
  semanaSanta <- sort(c(juevesSanto, viernesSanto, sabadoSanto, domingoResurrecion, lunesPascua))
  
  # Se pone en formato data.frame y se añade un indicador
  semanaSanta <- data.frame(fechas=semanaSanta, semanaSanta=rep(1,length(semanaSanta)))
  
  # Se añaden a la tabla maestra de fechas
  todasLasFechas_2 <- merge(x = todasLasFechas, y = semanaSanta, by = "fechas", all.x = TRUE)
  
  # Se reemplazan los NAs por 0
  todasLasFechas_2$semanaSanta[is.na(todasLasFechas_2$semanaSanta)] <- 0
  
  
  ######################################
  #     Cálculo de la variable dt      #
  ######################################
  
  # 1. Definición de festivos:
  ############################
  
  calendario <- todasLasFechas
  
  calendario$diaSemana <- as.factor(wday(calendario$fecha))
  calendario$diaMes <- as.factor(day(calendario$fecha))
  calendario$mes <- as.factor(month(calendario$fecha))
  calendario$anyo <- as.factor(year(calendario$fecha))
  
  calendario$p_01ene <- ifelse(calendario$diaMes==1 & calendario$mes==1, 1, 0)
  calendario$p_06ene <- ifelse(calendario$diaMes==6 & calendario$mes==1, 1, 0)
  calendario$p_19mar <- ifelse(calendario$diaMes==19 & calendario$mes==3, 1, 0)
  calendario$p_01may <- ifelse(calendario$diaMes==1 & calendario$mes==5, 1, 0)
  calendario$p_15ago <- ifelse(calendario$diaMes==15 & calendario$mes==8, 1, 0)
  calendario$p_12oct <- ifelse(calendario$diaMes==12 & calendario$mes==10,1, 0)
  calendario$p_01nov <- ifelse(calendario$diaMes==1 & calendario$mes==11, 1 ,0)
  calendario$p_06dic <- ifelse(calendario$diaMes==6 & calendario$mes==12, 1 ,0)
  calendario$p_08dic <- ifelse(calendario$diaMes==8 & calendario$mes==12, 1 ,0)
  calendario$p_25dic <- ifelse(calendario$diaMes==25 & calendario$mes==12, 1 ,0)
  
  # Fiestas regionales: día de Extremadura
  calendario$p_08sep <- ifelse(calendario$diaMes==8 & calendario$mes==9, 1 ,0)
  # Fiestas locales: San Jorge
  calendario$p_23abr <- ifelse(calendario$diaMes==23 & calendario$mes==4, 1 ,0)
  
  calendario$festivo <- rowSums(subset(calendario, select=p_01ene:p_25dic))
  
  # La definición de la variable dt varía según la opción domingoYFestivosJuntos.
  
  if (domingoYFestivosJuntos==0){
    
    calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
    calendario$domingo <- ifelse(calendario$diaSemana==1, 1 ,0)
    
    # Días laborables: todos menos sábados y domingos
    calendario$laborable <- 1-calendario$sabado-calendario$domingo
    
  } else {
    
    calendario$sabado <- ifelse(calendario$diaSemana==7, 1 ,0)
    # Domingo=1 si domingo=1 o festivo=1
    calendario$domingo <- ifelse(calendario$domingo==1 | calendario$festivo==1, 1 ,0)
    
    # Días laborables: todos menos sábados y domingos/festivos
    calendario$laborable <- 1-calendario$sabado-calendario$domingo    
  }
  
  
  # 2. Definición de variable dt:
  ###############################
  
  # Se filtran las columnas de inter?s y se añade la Semana Santa
  
  calendario_2 <- calendario[, c("fechas", "mes", "anyo", "sabado", "domingo", "laborable", "festivo")]
  
  todasLasFechasFinal <- merge(x = todasLasFechas_2, y = calendario_2,
                               by = "fechas", all.x = TRUE)
  
  # Agregamos la serie a nivel a?o-mes
  
  calendarioAnyoMes <- aggregate(todasLasFechasFinal[,c("sabado","domingo",
                                                        "laborable", "semanaSanta", "festivo")],
                                 by=list(mes=todasLasFechasFinal$mes,
                                         anyo=todasLasFechasFinal$anyo),
                                 "sum")
  
  # Se calcula la variable dt:
  
  calendarioAnyoMes$dt <- calendarioAnyoMes$laborable-(5/2)*(calendarioAnyoMes$sabado+calendarioAnyoMes$domingo)
  
  
  ######################################
  #     Cálculo de años bisiestos      #
  ######################################
  
  calendarioAnyoMes$anyoNum <- as.numeric(levels(calendarioAnyoMes$anyo))[calendarioAnyoMes$anyo]
  
  calendarioAnyoMes$bisiesto <- ifelse(calendarioAnyoMes$mes==2 &(calendarioAnyoMes$anyoNum %% 4)==0, 1 ,0)
  
  
  #######################################################
  #     Tabla final con explicativas de calendario      #
  #######################################################
  
  if (domingoYFestivosJuntos==0){
    explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto", "festivo")])
  } else {
    explicativasCalendario <- cbind(fecha=variableFecha, calendarioAnyoMes[, c("semanaSanta", "dt", "bisiesto")])
  }
  
  return(explicativasCalendario)
  
}

Llamamos a la función que hemos creado con la opción de que domingos y días festivos se cuenten por separado.

explicativasCalendarioTrain <- calculoExplicativasCalendarioCaceres(datos.train$FECHA, domingoYFestivosJuntos=0)

calendarioTrain <- 
  as.matrix(
    explicativasCalendarioTrain[,c("semanaSanta", "dt", "bisiesto", "festivo")]
  )

Añadimos al ajuste las nuevas variables explicativas.

ajuste4conCalendario <- Arima(datos.train.ts,
                              order = c(0,1,1),
                              seasonal = list(order = c(1,0,0), period = 12),
                              xreg = calendarioTrain,
                              method = "ML")
ajuste4conCalendario
## Series: datos.train.ts 
## Regression with ARIMA(0,1,1)(1,0,0)[12] errors 
## 
## Coefficients:
##           ma1    sar1  semanaSanta      dt  bisiesto  festivo
##       -0.7200  0.1969      -0.0116  0.0150   -0.0227  -0.0614
## s.e.   0.0483  0.0734       0.0118  0.0056    0.1049   0.0212
## 
## sigma^2 = 0.0536:  log likelihood = 10.88
## AIC=-7.76   AICc=-7.15   BIC=15.01
coeftest(ajuste4conCalendario)
## 
## z test of coefficients:
## 
##               Estimate Std. Error  z value  Pr(>|z|)    
## ma1         -0.7199518  0.0483371 -14.8944 < 2.2e-16 ***
## sar1         0.1968911  0.0733748   2.6834  0.007289 ** 
## semanaSanta -0.0115610  0.0118108  -0.9788  0.327654    
## dt           0.0150418  0.0055963   2.6878  0.007192 ** 
## bisiesto    -0.0227074  0.1049431  -0.2164  0.828693    
## festivo     -0.0614426  0.0211852  -2.9003  0.003729 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vemos que el año bisiesto tiene un p-valor bastante grande.Además, el valor que puede aportarnos tener un día más en febrero ya está contemplado con la variable dt. Por otro lado, podemos mantener la Semana Santa ya que le podemos ver sentido de negocio.

calendarioTrain <- 
  as.matrix(
    explicativasCalendarioTrain[,c("semanaSanta", "dt", "festivo")]
  )

ajuste5conCalendario <- Arima(datos.train.ts,
                              order = c(0,1,1),
                              seasonal = list(order = c(1,0,0), period = 12),
                              xreg = calendarioTrain,
                              method = "ML")
ajuste5conCalendario
## Series: datos.train.ts 
## Regression with ARIMA(0,1,1)(1,0,0)[12] errors 
## 
## Coefficients:
##           ma1    sar1  semanaSanta      dt  festivo
##       -0.7203  0.1976      -0.0114  0.0150  -0.0608
## s.e.   0.0482  0.0733       0.0118  0.0056   0.0210
## 
## sigma^2 = 0.05333:  log likelihood = 10.86
## AIC=-9.71   AICc=-9.26   BIC=9.8
coeftest(ajuste5conCalendario)
## 
## z test of coefficients:
## 
##               Estimate Std. Error  z value  Pr(>|z|)    
## ma1         -0.7202627  0.0482284 -14.9344 < 2.2e-16 ***
## sar1         0.1975910  0.0732627   2.6970  0.006996 ** 
## semanaSanta -0.0113950  0.0117906  -0.9664  0.333821    
## dt           0.0150369  0.0055984   2.6859  0.007233 ** 
## festivo     -0.0607887  0.0209877  -2.8964  0.003775 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vemos correlaciones y vemos que ninguna variable parece tener correlación:

cor.arma(ajuste5conCalendario)
##                      ma1          sar1 semanaSanta            dt     festivo
## ma1          1.000000000 -0.1395694285 -0.02512203 -0.0060496156 -0.07048090
## sar1        -0.139569429  1.0000000000 -0.07245826 -0.0004302792 -0.03481278
## semanaSanta -0.025122031 -0.0724582600  1.00000000  0.0491736995  0.16027535
## dt          -0.006049616 -0.0004302792  0.04917370  1.0000000000  0.00428523
## festivo     -0.070480901 -0.0348127821  0.16027535  0.0042852305  1.00000000

Y el text de Ljung-Box para ruido blanco.

Box.test.2(residuals(ajuste5conCalendario),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard   p-value
## [1,]      6 0.7984990
## [2,]     12 0.9116931
## [3,]     18 0.8911888
## [4,]     24 0.7504270
## [5,]     30 0.9289244
## [6,]     36 0.8699718
## [7,]     42 0.8579565
## [8,]     48 0.9146346

Como resultado obtenemos que tenemos ruido blanco.

Identificación de outliers

Pasamos a identificar outliers con la función locate.ouliers.

Identificamos outliers de los tipos:

listaOutliersTrain <- locate.outliers(ajuste5conCalendario$residuals,
                                      pars = coefs2poly(ajuste5conCalendario),
                                      types = c("AO", "LS", "TC"),
                                      cval=3)

listaOutliersTrain$abststat=abs(listaOutliersTrain$tstat)

# Cruzamos con la tabla original para obtener la fecha

datos.train$ind <- as.numeric(rownames(datos.train))
listaOutliersTrainFecha <- merge(listaOutliersTrain, datos.train[,c("ind", "FECHA")], by = "ind")

arrange(listaOutliersTrainFecha,desc(listaOutliersTrainFecha$abststat))
##   ind type    coefhat     tstat abststat      FECHA
## 1  99   LS -0.4456849 -3.135621 3.135621 2011-03-01
## 2 126   TC -0.5369467 -3.119897 3.119897 2013-06-01
## 3  88   TC  0.5255352  3.053591 3.053591 2010-04-01

Nos han salido 2 outliers de tipo TC y 1 de tipo LS. Estos outliers los añadimos a las variables explicativas que ya teníamos.

outliersTrain <- outliers(c( "TC","LS","TC"), c(88, 99, 126))
outliersVariablesTrain <- outliers.effects(outliersTrain, length(ajuste5conCalendario$residuals))
calendarioMasOutliersTrain <- as.matrix(cbind(calendarioTrain, outliersVariablesTrain))

Ajustamos un nuevo modelo con todo ello.

ajuste5conCalendarioYOutliers <- Arima(datos.train.ts,
                                               order = c(0,1,1),
                                               seasonal = list(order = c(1,0,0), period = 12),
                                               xreg = calendarioMasOutliersTrain,
                                               method = "ML")

coeftest(ajuste5conCalendarioYOutliers)
## 
## z test of coefficients:
## 
##               Estimate Std. Error  z value  Pr(>|z|)    
## ma1         -0.7938138  0.0404222 -19.6381 < 2.2e-16 ***
## sar1         0.1901647  0.0759103   2.5051 0.0122409 *  
## semanaSanta -0.0141030  0.0114148  -1.2355 0.2166407    
## dt           0.0139374  0.0054566   2.5542 0.0106425 *  
## festivo     -0.0633449  0.0197687  -3.2043 0.0013539 ** 
## TC88         0.4022273  0.1792464   2.2440 0.0248330 *  
## LS99        -0.3986618  0.1380351  -2.8881 0.0038755 ** 
## TC126       -0.5954287  0.1754221  -3.3943 0.0006881 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si miramos las correlaciones vemos que no existen valores altos:

cor.arma(ajuste5conCalendarioYOutliers)
##                      ma1         sar1 semanaSanta           dt     festivo
## ma1          1.000000000 -0.102755031 -0.02714986  0.005189159 -0.01367233
## sar1        -0.102755031  1.000000000 -0.03635183 -0.004683842 -0.01923234
## semanaSanta -0.027149858 -0.036351834  1.00000000  0.053407042  0.15537449
## dt           0.005189159 -0.004683842  0.05340704  1.000000000  0.00851406
## festivo     -0.013672329 -0.019232340  0.15537449  0.008514060  1.00000000
## TC88         0.139030350 -0.008109016 -0.12474959 -0.013125127  0.05765664
## LS99         0.149116326  0.182599302 -0.07731113 -0.025581756  0.04218715
## TC126        0.164535310  0.134369279  0.03476912  0.074418665  0.06490375
##                     TC88        LS99      TC126
## ma1          0.139030350  0.14911633 0.16453531
## sar1        -0.008109016  0.18259930 0.13436928
## semanaSanta -0.124749593 -0.07731113 0.03476912
## dt          -0.013125127 -0.02558176 0.07441866
## festivo      0.057656643  0.04218715 0.06490375
## TC88         1.000000000  0.27256859 0.02299857
## LS99         0.272568590  1.00000000 0.05621866
## TC126        0.022998565  0.05621866 1.00000000

Volvemos a realizar el test de Ljung-Box y vemos que nos vuelve a dar ruido blanco.

Box.test.2(residuals(ajuste5conCalendarioYOutliers),
           nlag = c(6,12,18,24,30,36,42,48),
           type="Ljung-Box")
##      Retard   p-value
## [1,]      6 0.8503906
## [2,]     12 0.8326711
## [3,]     18 0.8147331
## [4,]     24 0.7642923
## [5,]     30 0.9193775
## [6,]     36 0.9386808
## [7,]     42 0.9294358
## [8,]     48 0.9033291

Predicción

Una vez ya obtenido nuestro modelo, vamos a pasar a intentar predecir el número de hipotecas que se concedieron en el año 2019.

Lo primero será añadir las variables explicativas del calendario a test.

explicativasCalendarioTest <- calculoExplicativasCalendarioCaceres(datos.test$FECHA,domingoYFestivosJuntos=0)
calendarioTest <- as.matrix(explicativasCalendarioTest[,c("semanaSanta", "dt", "festivo")])

Para los outliers tomaremos los del último año y los calcaremos en el año que vamos a predecir. Tras esto, juntamos las variables del calendario y las de los outliers.

outliersVariablesTest <- outliersVariablesTrain[(nrow(outliersVariablesTrain) - 11):nrow(outliersVariablesTrain),]
calendarioMasOutliersTest <- as.matrix(cbind(calendarioTest, outliersVariablesTest))
calendarioMasOutliersTest
##       semanaSanta   dt festivo         TC88 LS99        TC126
##  [1,]           0  3.0       2 3.927514e-15    1 3.022680e-09
##  [2,]           0  0.0       0 2.749260e-15    1 2.115876e-09
##  [3,]           0 -4.0       1 1.924482e-15    1 1.481113e-09
##  [4,]           5  2.0       0 1.347137e-15    1 1.036779e-09
##  [5,]           0  3.0       1 9.429961e-16    1 7.257455e-10
##  [6,]           0 -5.0       0 6.600972e-16    1 5.080219e-10
##  [7,]           0  3.0       0 4.620681e-16    1 3.556153e-10
##  [8,]           0 -0.5       1 3.234477e-16    1 2.489307e-10
##  [9,]           0 -1.5       0 2.264134e-16    1 1.742515e-10
## [10,]           0  3.0       1 1.584893e-16    1 1.219760e-10
## [11,]           0 -1.5       1 1.109425e-16    1 8.538323e-11
## [12,]           0 -0.5       3 7.765978e-17    1 5.976826e-11

Procedemos por lo tanto, a hacer la predicción de este nuevo año. Además, marcaremos los límites de confianza al 95%.

prediccion <- as.data.frame(predict(ajuste5conCalendarioYOutliers, 
                                    n.ahead=12,
                                    newxreg=calendarioMasOutliersTest))

#Límites de confianza

U <- exp(prediccion$pred + 2*prediccion$se)
L <- exp(prediccion$pred - 2*prediccion$se)

datos.pred <- data.frame(FECHA = datos.test$FECHA, 
                         Prediccion = exp(prediccion$pred + 0.5 * prediccion$se^2),
                         LimSup = U, 
                         LimInf = L)

Para ver nuestro resultado, elaboramos una gráfica mostrando la predicción, los límites de confianza y el valor real de la serie para el año 2019.

datos.real.pred <- merge(datos, datos.pred, by = "FECHA", all.x = T)

datos.real.pred$REAL <- datos.real.pred$TARGET

datos.real.pred$TARGET[datos.real.pred$FECHA>as.Date('01/12/2018',format='%d/%m/%Y')] <- NA

grafico <- ggplot(data = datos.real.pred) +
  geom_line(aes(x= FECHA, y = REAL), color = 'blue',
            alpha = 0.4, size = 0.8) +
  geom_line(aes(x= FECHA, y = TARGET), color = 'blue',
            alpha = 0.8, size = 0.8) +
  geom_line(aes(x= FECHA, y = Prediccion), color = 'darkred',
            size = 1)   +
  geom_line(aes(x= FECHA, y = LimSup), color = 'orange',
            size = 1)  +
  geom_line(aes(x= FECHA, y = LimInf), color = 'orange',
            size = 1) +
  xlab('FECHA') + ylab('Matriculaciones')

ggplotly(grafico)

Ajuste Automático

Podemos, además, comparar nuestro modelo con el ajuste que se obtiene al hacer un auto arima. Veamos el resultado.

ajusteAutomatico <- auto.arima(datos.train.ts,
                               max.d=1, max.D=1,
                               max.p=2, max.P=2,
                               max.q=2, max.Q=2, 
                               seasonal=TRUE,
                               ic="aic",
                               allowdrift=FALSE,
                               xreg=calendarioMasOutliersTrain,
                               stepwise=TRUE)
ajusteAutomatico
## Series: datos.train.ts 
## Regression with ARIMA(2,0,0) errors 
## 
## Coefficients:
##          ar1     ar2  intercept  semanaSanta      dt  festivo    TC88     LS99
##       0.3359  0.2129     6.1811      -0.0118  0.0144  -0.0766  0.3343  -1.2018
## s.e.  0.0730  0.0744     0.0558       0.0116  0.0051   0.0186  0.2375   0.0765
##         TC126
##       -0.6586
## s.e.   0.2284
## 
## sigma^2 = 0.05661:  log likelihood = 7.7
## AIC=4.6   AICc=5.82   BIC=37.18

El ajuste automático resuelve esta serie temporal con un ARIMA(2,0,0).

Vamos a calcular sus predicciones y a compararlas con las nuestras.

prediccion.auto <- as.data.frame(predict(ajusteAutomatico, 
                                    n.ahead=12,
                                    newxreg=calendarioMasOutliersTest))

datos.pred.auto <- data.frame(FECHA = datos.test$FECHA, 
                         Prediccion.auto = exp(prediccion.auto$pred + 0.5 * prediccion.auto$se^2))

datos.real.pred.auto <- merge(datos.real.pred, datos.pred.auto, by = "FECHA", all.x = T)

grafico <- ggplot(data = datos.real.pred.auto) +
  geom_line(aes(x= FECHA, y = REAL), color = 'blue',
            alpha = 0.4, size = 0.8) +
  geom_line(aes(x= FECHA, y = TARGET), color = 'blue',
            alpha = 0.8, size = 0.8) +
  geom_line(aes(x= FECHA, y = Prediccion), color = 'darkred',
            size = 1)   +
  geom_line(aes(x= FECHA, y = Prediccion.auto), color = 'orange',
            size = 1) +
  xlab('FECHA') + ylab('Matriculaciones') + ggtitle('Ajuste Manual vs Ajuste Automático')

ggplotly(grafico)

Vemos que nos arrojan unos resultados parecidos.